#-- SET ME   Tutorial should work with most binary black hole events
eventname = 'GW150914' 

# want plots?
make_plots = 1
plottype = "png"

# Standard python numerical analysis imports:
import numpy as np
import scipy
from scipy import signal
from scipy.interpolate import interp1d
from scipy.signal import butter, filtfilt, iirdesign, zpk2tf, freqz
import h5py
import json

# the IPython magic below must be commented out in the .py file, since it doesn't work there.
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

# LIGO-specific readligo.py 
import LOSC_Event_tutorial.readligo as rl

from pprint import pprint

# you might get a matplotlib warning here; you can ignore it.

# Read the event properties from a local json file
fnjson = "LOSC_Event_tutorial/BBH_events_v3.json"
try:
    events = json.load(open(fnjson,"r"))
except IOError:
    print("Cannot find resource file "+fnjson)
    print("You can download it from https://losc.ligo.org/s/events/"+fnjson)
    print("Quitting.")
    quit()

# did the user select the eventname ?
try: 
    events[eventname]
except:
    print('You must select an eventname that is in '+fnjson+'! Quitting.')
    quit()

# Extract the parameters for the desired event:
event = events[eventname]
fn_H1 = event['fn_H1']              # File name for H1 data
fn_L1 = event['fn_L1']              # File name for L1 data
fn_template = event['fn_template']  # File name for template waveform
fs = event['fs']                    # Set sampling rate
tevent = event['tevent']            # Set approximate event GPS time
fband = event['fband']              # frequency band for bandpassing signal
print("Reading in parameters for event " + event["name"])
pprint(event)

Reading in parameters for event GW150914
{'a1': 0.355,
 'a2': -0.769,
 'approx': 'lalsim.SEOBNRv2',
 'f_min': 10.0,
 'fband': [43.0, 300.0],
 'fn_H1': 'H-H1_LOSC_4_V2-1126259446-32.hdf5',
 'fn_L1': 'L-L1_LOSC_4_V2-1126259446-32.hdf5',
 'fn_template': 'GW150914_4_template.hdf5',
 'fs': 4096,
 'm1': 41.743,
 'm2': 29.237,
 'name': 'GW150914',
 'tevent': 1126259462.44,
 'utcevent': '2015-09-14T09:50:45.44'}

# 确认要加载的数据的存在性
!ls -lht LOSC_Event_tutorial/$fn_H1
!ls -lht LOSC_Event_tutorial/$fn_L1

-rw-r--r-- 1 root root 1017K Dec 12 12:35 LOSC_Event_tutorial/H-H1_LOSC_4_V2-1126259446-32.hdf5
-rw-r--r-- 1 root root 984K Dec 12 12:35 LOSC_Event_tutorial/L-L1_LOSC_4_V2-1126259446-32.hdf5

#----------------------------------------------------------------
# Load LIGO data from a single file.
# FIRST, define the filenames fn_H1 and fn_L1, above.
#----------------------------------------------------------------
try:
    # read in data from H1 and L1, if available:
    strain_H1, time_H1, chan_dict_H1 = rl.loaddata('LOSC_Event_tutorial/'+fn_H1, 'H1')
    strain_L1, time_L1, chan_dict_L1 = rl.loaddata('LOSC_Event_tutorial/'+fn_L1, 'L1')
except:
    print("Cannot find data files!")
    print("You can download them from https://losc.ligo.org/s/events/"+eventname)
    print("Quitting.")
    # quit()  # kernel restart!

# 了解你的数据！
print(strain_H1.shape, time_H1.shape, strain_H1.dtype, time_H1.dtype)
print(strain_L1.shape, time_L1.shape, strain_L1.dtype, time_L1.dtype)

(131072,) (131072,) float64 float64
(131072,) (131072,) float64 float64

# 数据一致性检查
assert np.alltrue(time_H1 == time_L1)
assert np.all(np.diff(time_H1))

# 采样率
print(1/np.diff(time_H1)[0], 'Hz')

# 采样频率（sample frequency）物理意义：物理上的每一秒时间范围里，有多少个采样点。

4096.0 Hz

# 测试一下。。。
time_H1[4096] - time_H1[0]

1.0

# pprint(chan_dict_H1)
pprint(chan_dict_L1)
chan_dict_H1['BURST_CAT1'].shape

{'BURST_CAT1': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32),
 'BURST_CAT2': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32),
 'BURST_CAT3': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32),
 'CBC_CAT1': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32),
 'CBC_CAT2': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32),
 'CBC_CAT3': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32),
 'DATA': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32),
 'DEFAULT': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32),
 'NO_BURST_HW_INJ': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32),
 'NO_CBC_HW_INJ': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32),
 'NO_CW_HW_INJ': array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=uint32),
 'NO_DETCHAR_HW_INJ': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32),
 'NO_STOCH_HW_INJ': array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint32)}

(32,)

# both H1 and L1 will have the same time vector, so:
time = time_H1
# the time sample interval (uniformly sampled!)
dt = time[1] - time[0]

# Let's look at the data and print out some stuff:

print('time_H1: len, min, mean, max = ', \
    len(time_H1), time_H1.min(), time_H1.mean(), time_H1.max() )
print('strain_H1: len, min, mean, max = ', \
    len(strain_H1), strain_H1.min(),strain_H1.mean(),strain_H1.max())
print( 'strain_L1: len, min, mean, max = ', \
    len(strain_L1), strain_L1.min(),strain_L1.mean(),strain_L1.max())

#What's in chan_dict?  (See also https://losc.ligo.org/tutorials/)
bits = chan_dict_H1['DATA']
print("For H1, {0} out of {1} seconds contain usable DATA".format(bits.sum(), len(bits)))
bits = chan_dict_L1['DATA']
print("For L1, {0} out of {1} seconds contain usable DATA".format(bits.sum(), len(bits)))

time_H1: len, min, mean, max =  131072 1126259446.0 1126259461.999878 1126259477.9997559
strain_H1: len, min, mean, max =  131072 -7.044665943156067e-19 5.895522509246437e-23 7.706262192397465e-19
strain_L1: len, min, mean, max =  131072 -1.8697138664279764e-18 -1.0522332249909908e-18 -4.60035111311666e-20
For H1, 32 out of 32 seconds contain usable DATA
For L1, 32 out of 32 seconds contain usable DATA

# plot +- deltat seconds around the event:
# index into the strain time series for this time interval:
deltat = 5
indxt = np.where((time >= tevent-deltat) & (time < tevent+deltat))  # 切片 GW150914 附近 ± 5秒的索引
print('The event at', tevent)

if make_plots:
    plt.figure()
    plt.plot(time[indxt]-tevent, strain_H1[indxt],'r',label='H1 strain')
    plt.plot(time[indxt]-tevent, strain_L1[indxt],'g',label='L1 strain')
    plt.axvline(0, color='k', ls='--')
    plt.xlabel('time (s) since '+str(tevent))
    plt.ylabel('strain')
    plt.legend(loc='lower right')
    plt.title('Advanced LIGO strain data near '+eventname)
    # plt.savefig(eventname+'_strain.'+plottype)

The event at 1126259462.44

import numpy as np
def DFT_slow(x):
    """Compute the discrete Fourier Transform of the 1D array x"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)

fs = 4096
time_array = np.arange(0, 1, 1/fs)
x = np.sin(time_array)
x

array([0.00000000e+00, 2.44140623e-04, 4.88281231e-04, ...,
       8.41075030e-01, 8.41207065e-01, 8.41339050e-01])

N = x.shape[0]
N

4096

n = np.arange(N)
k = n.reshape((N, 1))

n.shape, k.shape, (k * n).shape

((4096,), (4096, 1), (4096, 4096))

(k * n)

array([[       0,        0,        0, ...,        0,        0,        0],
       [       0,        1,        2, ...,     4093,     4094,     4095],
       [       0,        2,        4, ...,     8186,     8188,     8190],
       ...,
       [       0,     4093,     8186, ..., 16752649, 16756742, 16760835],
       [       0,     4094,     8188, ..., 16756742, 16760836, 16764930],
       [       0,     4095,     8190, ..., 16760835, 16764930, 16769025]])

type(1j)

complex

M = np.exp(-2j * np.pi * k * n / N)
print(M.shape)
M

(4096, 4096)

array([[1.        +0.j        , 1.        +0.j        ,
        1.        +0.j        , ..., 1.        +0.j        ,
        1.        +0.j        , 1.        +0.j        ],
       [1.        +0.j        , 0.99999882-0.00153398j,
        0.99999529-0.00306796j, ..., 0.99998941+0.00460193j,
        0.99999529+0.00306796j, 0.99999882+0.00153398j],
       [1.        +0.j        , 0.99999529-0.00306796j,
        0.99998118-0.00613588j, ..., 0.99995764+0.00920375j,
        0.99998118+0.00613588j, 0.99999529+0.00306796j],
       ...,
       [1.        +0.j        , 0.99998941+0.00460193j,
        0.99995764+0.00920375j, ..., 0.9999047 -0.01380539j,
        0.99995764-0.00920375j, 0.99998941-0.00460193j],
       [1.        +0.j        , 0.99999529+0.00306796j,
        0.99998118+0.00613588j, ..., 0.99995764-0.00920375j,
        0.99998118-0.00613588j, 0.99999529-0.00306796j],
       [1.        +0.j        , 0.99999882+0.00153398j,
        0.99999529+0.00306796j, ..., 0.99998941-0.00460193j,
        0.99999529-0.00306796j, 0.99999882-0.00153398j]])

np.dot(M, x).shape

(4096,)

def DFT_slow(x):
    """Compute the discrete Fourier Transform of the 1D array x"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)

DFT_slow(x)

array([1882.50101032  +0.j        ,  -49.3552352 +562.80983113j,
        -12.42047533+276.02463083j, ...,   -5.73514268-183.36702601j,
        -12.42047533-276.02463083j,  -49.3552352 -562.80983113j])

# Compute the one-dimensional discrete Fourier Transform.

# This function computes the one-dimensional *n*-point discrete Fourier
# Transform (DFT) with the efficient Fast Fourier Transform (FFT)
# algorithm [CT].

np.fft.fft(x)

array([1882.50101032  +0.j        ,  -49.3552352 +562.80983113j,
        -12.42047533+276.02463083j, ...,   -5.73514268-183.36702601j,
        -12.42047533-276.02463083j,  -49.3552352 -562.80983113j])

print(np.allclose(DFT_slow(x), np.fft.fft(x)))

True

def FFT(x):
    """A recursive implementation of the 1D Cooley-Tukey FFT"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    
    if N % 2 > 0:
        raise ValueError("size of x must be a power of 2")
    elif N <= 32:  # this cutoff should be optimized
        return DFT_slow(x)
    else:
        X_even = FFT(x[::2])
        X_odd = FFT(x[1::2])
        factor = np.exp(-2j * np.pi * np.arange(N) / N)
        return np.concatenate([X_even + factor[:N // 2] * X_odd,
                               X_even + factor[N // 2:] * X_odd])

def FFT(x):
    """A recursive implementation of the 1D Cooley-Tukey FFT"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    
    if N % 2 > 0:
        raise ValueError("size of x must be a power of 2")
    elif N <= 32:  # this cutoff should be optimized
        return DFT_slow(x)
    else:
        X_even = FFT(x[::2])
        X_odd = FFT(x[1::2])
        factor = np.exp(-2j * np.pi * np.arange(N) / N)
        return np.concatenate([X_even + factor[:N // 2] * X_odd,
                               X_even + factor[N // 2:] * X_odd])

print(np.allclose(DFT_slow(x), FFT(x)))

True

def FFT_vectorized(x):
    """A vectorized, non-recursive version of the Cooley-Tukey FFT"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]

    if np.log2(N) % 1 > 0:
        raise ValueError("size of x must be a power of 2")

    # N_min here is equivalent to the stopping condition above,
    # and should be a power of 2
    N_min = min(N, 32)
    
    # Perform an O[N^2] DFT on all length-N_min sub-problems at once
    n = np.arange(N_min)
    k = n[:, None]
    M = np.exp(-2j * np.pi * n * k / N_min)
    X = np.dot(M, x.reshape((N_min, -1)))

    # build-up each level of the recursive calculation all at once
    while X.shape[0] < N:
        X_even = X[:, :X.shape[1] // 2]
        X_odd = X[:, X.shape[1] // 2:]
        factor = np.exp(-1j * np.pi * np.arange(X.shape[0])
                        / X.shape[0])[:, None]
        X = np.vstack([X_even + factor * X_odd,
                       X_even - factor * X_odd])

    return X.ravel()

def FFT_vectorized(x):
    """A vectorized, non-recursive version of the Cooley-Tukey FFT"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]

    if np.log2(N) % 1 > 0:
        raise ValueError("size of x must be a power of 2")

    # N_min here is equivalent to the stopping condition above,
    # and should be a power of 2
    N_min = min(N, 32)
    
    # Perform an O[N^2] DFT on all length-N_min sub-problems at once
    n = np.arange(N_min)
    k = n[:, None]
    M = np.exp(-2j * np.pi * n * k / N_min)
    X = np.dot(M, x.reshape((N_min, -1)))

    # build-up each level of the recursive calculation all at once
    while X.shape[0] < N:
        X_even = X[:, :X.shape[1] // 2]
        X_odd = X[:, X.shape[1] // 2:]
        factor = np.exp(-1j * np.pi * np.arange(X.shape[0])
                        / X.shape[0])[:, None]
        X = np.vstack([X_even + factor * X_odd,
                       X_even - factor * X_odd])

    return X.ravel()

print(np.allclose(DFT_slow(x), FFT_vectorized(x)))

True

x = np.random.random(1024)
print(np.allclose(DFT_slow(x), np.fft.fft(x)))
print(np.allclose(DFT_slow(x), FFT(x)))
print(np.allclose(DFT_slow(x), FFT_vectorized(x)))

import scipy
print(np.allclose(DFT_slow(x), scipy.fftpack.fft(x)))

True
True
True
True

x = np.random.random(1024)
%timeit DFT_slow(x)
%timeit FFT(x)
%timeit FFT_vectorized(x)
%timeit np.fft.fft(x)
%timeit scipy.fftpack.fft(x)

59 ms ± 4.79 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
2.41 ms ± 43.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
261 µs ± 339 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
14 µs ± 15.2 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
12.3 µs ± 19.7 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

# Signature: np.fft.fftfreq(n, d=1.0)
# Docstring:
# Return the Discrete Fourier Transform sample frequencies.

# The returned float array `f` contains the frequency bin centers in cycles
# per unit of the sample spacing (with zero at the start).  For instance, if
# the sample spacing is in seconds, then the frequency unit is cycles/second.

# Given a window length `n` and a sample spacing `d`::

#   f = [0, 1, ...,   n/2-1,     -n/2, ..., -1] / (d*n)   if n is even
#   f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n)   if n is odd

np.fft.fftfreq

N = 4096 # number of samples

fs = 100 # sampling frequency
dt = 1./fs # time discretisation

time_array = np.arange(N)*dt # time sampling

A = 0.067 # sine wave amplitude
f = 10.24 # sine wave frequency
strain = A*np.sin(2*np.pi*f*time_array) # sine wave

plt.plot(strain[:300])

[<matplotlib.lines.Line2D at 0x7fdcd7f30a30>]

freqs = np.fft.fftfreq(N, d=1./fs)
amp_spectrum_twosided = abs(DFT_slow(strain)/N)   # 用总采样的数目 N 来归一化 => Amplitude spectrum 幅度谱

# np.fft.fftshift: shift the zero-frequency component to the center of the spectrum.
assert np.alltrue(np.fft.fftshift(freqs) == np.roll(freqs, shift=N//2))

plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(np.fft.fftshift(freqs),
         np.fft.fftshift(amp_spectrum_twosided), label='Sine wave')
plt.xscale('symlog')
plt.yscale('log')
plt.axvline(f, label=f'{f}Hz', ls=':', color='k')
plt.title('Two-sided amplitude spectrum')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [$V_{RMS}$]')
plt.legend()
plt.grid();

plt.subplot(1,2,2)
plt.plot(time_array, strain, label='Sine wave')
plt.title('Time history')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude [V]')
plt.xlim(0,3)
plt.legend()
plt.grid();

fs

100

max(freqs)  # 果然是小于 fs=100Hz 的一半。

49.9755859375

def rms_flat(a):
    """
    Return the root mean square of all the elements of *a*, flattened out.
    Ref: https://matplotlib.org/2.2.5/_modules/matplotlib/mlab.html#rms_flat
    """
    return np.sqrt(np.mean(np.abs(a) ** 2))

print(f'{max(amp_spectrum_twosided):.6f}, {A/(2*np.sqrt(2)):.6f}')

0.024192, 0.023688

# 1. 截断FFT ： 基本上就是只计算感兴趣的项；
freqs_truncated = np.fft.rfftfreq(N, d=1./fs)
# Compute the one-dimensional discrete Fourier Transform for real input.
amp_spectrum_truncated = abs(np.fft.rfft(strain)/N)   # 用总采样的数目 N 来归一化

plt.figure(figsize=(6,4))
plt.loglog(freqs_truncated, amp_spectrum_truncated, label='Sine wave')
plt.axvline(f, label=f'{f}Hz', ls=':', color='k')
plt.title('Truncated amplitude spectrum')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [$V_{RMS}$]')
plt.legend()
plt.grid();

np.allclose(amp_spectrum_truncated, amp_spectrum_twosided[:amp_spectrum_truncated.shape[0]]  )

True

# 2. 折叠FFT ： 计算FFT的坐标，就好像它们围绕 Nyquist 频率折叠以获得 单边频谱（one-sided spectrum）。
scale = 2.*np.ones_like(freqs_truncated) # scale rfft components by a factor of 2
scale[0] = 1. #  the DC component is not scaled
if scale.size % 2 == True: # if there is a Nyquist component...
    scale[-1] = 1. # ...then it is not scaled
amp_spectrum_onesided = amp_spectrum_truncated * scale

plt.figure(figsize=(6,4))
plt.loglog(freqs_truncated, amp_spectrum_onesided, label='Sine wave')
plt.axvline(f, label=f'{f}Hz', ls=':', color='k')
plt.title('One-sided amplitude spectrum')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [$V_{RMS}$]')
plt.legend()
plt.grid();

print(f'{max(amp_spectrum_onesided):.6f}, {A/np.sqrt(2):.6f}')

0.048384, 0.047376

N = 4096 # number of samples

fs = 100 # sampling frequency
dt = 1./fs # time discretisation

time_array = np.arange(N)*dt # time sampling

A = 0.067 # sine wave amplitude
f = 10.24 # sine wave frequency
strain = A*np.sin(2*np.pi*f*time_array) # sine wave

freqs = np.fft.fftfreq(N, d=1./fs)
poser_spectrum_twosided = DFT_slow(strain) * np.conj(DFT_slow(strain))
poser_spectrum_twosided = np.real(poser_spectrum_twosided)

abs(DFT_slow(strain))**2 # 模平方

array([0.03146645, 0.03146687, 0.0314681 , ..., 0.03147016, 0.0314681 ,
       0.03146687])

plt.figure(figsize=(6,4))
plt.plot(np.fft.fftshift(freqs),
         np.fft.fftshift(poser_spectrum_twosided)/N**2, label='Sine wave')
plt.xscale('symlog')
plt.yscale('log')
plt.axvline(f, label=f'{f}Hz', ls=':', color='k')
plt.title('Two-sided power spectrum (fft)')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [$V_{RMS}^2$]')
plt.legend()
plt.grid();

print(f'{max(poser_spectrum_twosided)/N**2:.6f}, {(A**2/8):.6f}')

0.000585, 0.000561

# 紧跟上一小节的例子
win = np.hanning(strain.size)
plt.figure(figsize=(5,3))
plt.plot(win, 'r-', label='Hanning window')
plt.grid()
plt.legend()

freqs = np.fft.fftfreq(N, d=1./fs)
poser_spectrum_twosided_win = DFT_slow(strain * win) * np.conj(DFT_slow(strain * win))
poser_spectrum_twosided_win = np.real(poser_spectrum_twosided_win)

plt.figure(figsize=(6,4))
plt.plot(np.fft.fftshift(freqs),
         np.fft.fftshift(poser_spectrum_twosided)/N**2, label='FFT')
plt.plot(np.fft.fftshift(freqs),
         np.fft.fftshift(poser_spectrum_twosided_win)/N**2, label='FFT w. window')
plt.xscale('symlog')
plt.yscale('log')
plt.axvline(f, label=f'{f}Hz', ls=':', color='k')
plt.title('Two-sided power spectrum (fft)')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [$V_{RMS}^2$]')
plt.legend()
plt.grid();

print(f'{max(poser_spectrum_twosided_win)/(win**2).mean()/N**2:.6f}, {(A**2/8):.6f}')

0.000587, 0.000561

# First generate some sampled data
# define the sample rate
fs = 4096   # [Hz]
# generate a time vector
dt = 1./fs
ttot = 10.  # [sec]
t = np.arange(fs*ttot)*dt
# make a sine wave
f0 = 250.
signal = np.sin(2*np.pi*f0*t)
# make some random Gaussian noise
noise = np.random.randn(t.size)

# plot in time
plt.figure()
plt.plot(t,noise)
plt.xlabel('time (s)')
plt.ylabel('value (arb units)')
plt.grid()
plt.title('Random Gaussian data')

# plot the noise, sine wave, and their sum
strain = signal + noise
plt.figure()
plt.plot(t,strain,'r',label='sin+randn')
plt.plot(t,noise,'g', label='randn')
plt.plot(t,signal,'b',label='sin')
plt.xlabel('time (s)')
plt.ylabel('value (arb units)')
plt.grid()
plt.legend()

<matplotlib.legend.Legend at 0x7fdcd745aa10>

NFFT = fs//8    # 要解析的短时窗口序列的长度 (number of per segment, nperseg)
noverlap = NFFT//2  # 每个相邻短时窗口序列之间的重叠长度 (num of overlap)。一般来说通常取定为 nperseg//2，即 50% 的窗口重叠。

Pxx, freqs = mlab.psd(strain, Fs = fs, NFFT = NFFT, noverlap=NFFT//2)

plt.figure(figsize=(12, 4))
plt.subplot(1,2,1)
plt.plot(freqs, np.sqrt(Pxx))
plt.xlim([20.,fs/2.])
plt.grid()
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude Spectral Density (ASD)');

plt.subplot(1,2,2)
plt.plot(freqs, Pxx)
plt.xlim([20.,fs/2.])
plt.grid()
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (PSD)');

NFFT = nperseg = fs//8    # 要解析的短时窗口序列的长度 (number of per segment, nperseg)
noverlap = NFFT//2  # 每个相邻短时窗口序列之间的重叠长度 (num of overlap)。一般来说通常取定为 nperseg//2，即 50% 的窗口重叠。
nstep = NFFT - noverlap  # 间接计算出来的参数 (num of step)，为了之后计算窗口滑动数目带来便利。
win = np.hanning(NFFT)   # 对每个短时窗口序列加窗

# onesided
numFreqs = NFFT//2 + 1  # 由于 NFFT 是偶数
scaling_factor = 2.

NFFT, noverlap, nstep, win.size

(512, 256, 256, 512)

# 此时就可以把 freqs 坐标结果先算出来
freqs_scratch = np.fft.fftfreq(NFFT, 1/fs)[:numFreqs]

# get the last value correctly, it is negative otherwise
freqs_scratch[-1] *= -1
freqs_scratch.shape

(257,)

# Check
np.alltrue(freqs == freqs_scratch)

True

# 对 strain 数据进行切片
slices = [ (i*nstep , i*nstep+nperseg)
          for i in range(strain.size)
          if i*nstep+nperseg <= strain.size ]
result = np.asarray([strain[s[0]:s[1]] for s in slices])  # (159, 512) => (time, freq)

# Apply window by multiplication
result *= win  # (159, 512) x (512,)

# Perform the fft.
# 1-D *nfft*-point discrete Fourier Transform (DFT) of a real-valued array 
# result = scipy.fft.rfft(result.real, n=nfft)[:,:numFreqs]
result = np.fft.fft(result, n=NFFT)[:,:numFreqs]

result = np.conj(result) * result  # 取模平方 = | |^2

# Also include scaling factors for one-sided densities and dividing by
# the sampling frequency, if desired. Scale everything, except the DC
# component and the NFFT/2 component:
slc = slice(1, -1, None)
result[:,slc] *= scaling_factor

# The resulting density values should be scaled by the scaling
#     frequency, which gives density in units of 1/Hz.  This allows for
#     integration over the returned frequency values.  The default is True for
#     MATLAB compatibility.
result /= fs  # You know why...

# Energy correction.
# Scale the spectrum by the norm of the window to compensate for
# windowing loss; see Bendat & Piersol Sec 11.5.2.
result /= (win**2).sum()
result.shape

(159, 257)

# 手撸 one-sided PSD
Pxx_scratch = result.mean(axis=0).real
Pxx_scratch.shape

(257,)

# Check
np.alltrue(Pxx == Pxx_scratch)

True

# Check
plt.figure(figsize=(12, 4))
plt.subplot(1,2,1)
plt.plot(freqs_scratch, np.sqrt(Pxx_scratch))
plt.xlim([20.,fs/2.])
plt.grid()
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude Spectral Density (ASD, strain/rtHz)');

plt.subplot(1,2,2)
plt.plot(freqs_scratch, Pxx_scratch)
plt.xlim([20.,fs/2.])
plt.grid()
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (PSD)');

def spectral_helper(strain, fs, NFFT, noverlap):
    nperseg = NFFT
    nstep = NFFT - noverlap  # 间接计算出来的参数 (num of step)，为了之后计算窗口滑动数目带来便利。
    win = np.hanning(NFFT)   # 对每个短时窗口序列加窗

    # onesided
    numFreqs = NFFT//2 + 1  # 由于 NFFT 是偶数
    scaling_factor = 2.


    # onesided
    numFreqs = NFFT//2 + 1  # 由于 NFFT 是偶数
    scaling_factor = 2.


    # 对 strain 数据进行切片
    slices = [ (i*nstep , i*nstep+nperseg)
              for i in range(strain.size)
              if i*nstep+nperseg <= strain.size ]
    result = np.asarray([strain[s[0]:s[1]] for s in slices])  # (159, 512) => (time, freq)

    # Apply window by multiplication
    result *= win  # (159, 512) x (512,)

    # Perform the fft.
    # 1-D *nfft*-point discrete Fourier Transform (DFT) of a real-valued array 
    # result = scipy.fft.rfft(result.real, n=nfft)[:,:numFreqs]
    result = np.fft.fft(result, n=NFFT)[:,:numFreqs]

    result = np.conj(result) * result  # 取模平方 = | |^2

    # Also include scaling factors for one-sided densities and dividing by
    # the sampling frequency, if desired. Scale everything, except the DC
    # component and the NFFT/2 component:
    slc = slice(1, -1, None)
    result[:,slc] *= scaling_factor

    # The resulting density values should be scaled by the scaling
    #     frequency, which gives density in units of 1/Hz.  This allows for
    #     integration over the returned frequency values.  The default is True for
    #     MATLAB compatibility.
    result /= fs  # You know why...

    # Energy correction.
    # Scale the spectrum by the norm of the window to compensate for
    # windowing loss; see Bendat & Piersol Sec 11.5.2.
    result /= (win**2).sum()
    return result

result = spectral_helper(strain, fs, NFFT, noverlap)

# 手撸 one-sided PSD
Pxx_scratch = result.mean(axis=0).real

# Check
np.alltrue(Pxx == Pxx_scratch)

True

spec, freqs, bins, im = plt.specgram(strain, Fs = fs, NFFT = NFFT, noverlap=noverlap)
plt.xlabel('Time (s)')
plt.ylabel('Freq (Hz)')

spec.shape

(257, 159)

spec_scratch = spectral_helper(strain, fs, NFFT, noverlap).real
spec_scratch.shape

(159, 257)

# Check
np.alltrue(spec.T == spec_scratch)

True

make_psds = 1
if make_psds:
    # number of sample for the fast fourier transform:
    NFFT = 4*fs
    Pxx_H1, freqs = mlab.psd(strain_H1, Fs = fs, NFFT = NFFT)
    Pxx_L1, freqs = mlab.psd(strain_L1, Fs = fs, NFFT = NFFT)

    # We will use interpolations of the ASDs computed above for whitening:
    psd_H1 = interp1d(freqs, Pxx_H1)
    psd_L1 = interp1d(freqs, Pxx_L1)

    # Here is an approximate, smoothed PSD for H1 during O1, with no lines. We'll use it later.    
    Pxx = (1.e-22*(18./(0.1+freqs))**2)**2+0.7e-23**2+((freqs/2000.)*4.e-23)**2
    psd_smooth = interp1d(freqs, Pxx)

if make_plots:
    # plot the ASDs, with the template overlaid:
    f_min = 20.
    f_max = 2000. 
    plt.figure(figsize=(10,8))
    plt.loglog(freqs, np.sqrt(Pxx_L1),'g',label='L1 strain')
    plt.loglog(freqs, np.sqrt(Pxx_H1),'r',label='H1 strain')
    plt.loglog(freqs, np.sqrt(Pxx),'k',label='H1 strain, O1 smooth model')
    plt.axis([f_min, f_max, 1e-24, 1e-19])
    plt.grid('on')
    plt.ylabel('ASD (strain/rtHz)')
    plt.xlabel('Freq (Hz)')
    plt.legend(loc='upper center')
    plt.title('Advanced LIGO strain data near '+eventname)
    # plt.savefig(eventname+'_ASDs.'+plottype)

# function to whiten data
def whiten(strain, interp_psd, dt):
    Nt = len(strain)
    freqs = np.fft.rfftfreq(Nt, dt)
    freqs1 = np.linspace(0,2048,Nt//2+1)

    # whitening: transform to freq domain, divide by asd, then transform back, 
    # taking care to get normalization right.
    hf = np.fft.rfft(strain)
    norm = 1./np.sqrt(1./(dt*2))
    white_hf = hf / np.sqrt(interp_psd(freqs)) * norm
    white_ht = np.fft.irfft(white_hf, n=Nt)
    return white_ht

whiten_data = 1
if whiten_data:
    # now whiten the data from H1 and L1, and the template (use H1 PSD):
    strain_H1_whiten = whiten(strain_H1,psd_H1,dt)
    strain_L1_whiten = whiten(strain_L1,psd_L1,dt)
    
    # We need to suppress the high frequency noise (no signal!) with some bandpassing:
    bb, ab = butter(4, [fband[0]*2./fs, fband[1]*2./fs], btype='band')
    normalization = np.sqrt((fband[1]-fband[0])/(fs/2))
    strain_H1_whitenbp = filtfilt(bb, ab, strain_H1_whiten) / normalization
    strain_L1_whitenbp = filtfilt(bb, ab, strain_L1_whiten) / normalization

fband

[43.0, 300.0]

# plot +- deltat seconds around the event:
# index into the strain time series for this time interval:
deltat = 5
indxt = np.where((time >= tevent-deltat) & (time < tevent+deltat))  # 切片 GW150914 附近 ± 5秒的索引
print('The event at', tevent)

colors = ['#1f77b4', '#ff7f0e', '#2ca02c']

plt.figure(figsize=(12,12))
for idx, (strain_type, text) in enumerate(zip([strain_H1, strain_H1_whiten, strain_H1_whitenbp],
                               ['', ' (whiten)', ' (whiten+bandpass)'])):
    plt.subplot(3,2,(idx*2+1))
    plt.plot(time[indxt]-tevent, strain_type[indxt],f'{colors[idx]}',label='H1 strain'+text)
    plt.axvline(0, color='k', ls='--')
    plt.ylabel('strain')
    plt.legend(loc='lower right')
    if idx==0:
        plt.title('Advanced LIGO strain data near '+eventname)
    plt.xlabel('time (s) since '+str(tevent))

    plt.subplot(2,2,2)
    Pxx_type, freqs = mlab.psd(strain_type, Fs = fs, NFFT = NFFT)
    plt.loglog(freqs, np.sqrt(Pxx_type),f'{colors[idx]}',label='H1 strain'+text)
    # plt.axis([f_min, f_max, 1e-24, 1e-19])
    plt.grid('on')
    plt.ylabel('ASD (strain/rtHz)')
    plt.xlabel('Freq (Hz)')
    plt.legend(loc='upper center')
    plt.title('Advanced LIGO strain data near '+eventname)
    
plt.subplot(3,2,6)
plt.plot(time[indxt]-tevent, strain_type[indxt],f'{colors[idx]}',label='H1 strain'+text)
plt.axvline(0, color='k', ls='--')
plt.ylabel('strain')
plt.legend(loc='lower right')
plt.xlabel('time (s) since '+str(tevent))
plt.xlim(-0.5,0.1);

The event at 1126259462.44

if make_plots:
    # index into the strain time series for this time interval:
    indxt = np.where((time >= tevent-deltat) & (time < tevent+deltat))

    # pick a shorter FTT time interval, like 1/8 of a second:
    NFFT = int(fs/8)
    # and with a lot of overlap, to resolve short-time features:
    NOVL = int(NFFT*15./16)
    # and choose a window that minimizes "spectral leakage" 
    # (https://en.wikipedia.org/wiki/Spectral_leakage)
    window = np.blackman(NFFT)

    # the right colormap is all-important! See:
    # http://matplotlib.org/examples/color/colormaps_reference.html
    # viridis seems to be the best for our purposes, but it's new; if you don't have it, you can settle for ocean.
    #spec_cmap='viridis'
    spec_cmap='ocean'

    # Plot the H1 spectrogram:
    plt.figure(figsize=(10,6))
    spec_H1, freqs, bins, im = plt.specgram(strain_H1[indxt], NFFT=NFFT, Fs=fs, window=window, 
                                            noverlap=NOVL, cmap=spec_cmap, xextent=[-deltat,deltat])
    plt.xlabel('time (s) since '+str(tevent))
    plt.ylabel('Frequency (Hz)')
    plt.colorbar()
    plt.axis([-deltat, deltat, 0, 2000])
    plt.title('aLIGO H1 strain data near '+eventname)
    # plt.savefig(eventname+'_H1_spectrogram.'+plottype)

    # Plot the L1 spectrogram:
    plt.figure(figsize=(10,6))
    spec_H1, freqs, bins, im = plt.specgram(strain_L1[indxt], NFFT=NFFT, Fs=fs, window=window, 
                                            noverlap=NOVL, cmap=spec_cmap, xextent=[-deltat,deltat])
    plt.xlabel('time (s) since '+str(tevent))
    plt.ylabel('Frequency (Hz)')
    plt.colorbar()
    plt.axis([-deltat, deltat, 0, 2000])
    plt.title('aLIGO L1 strain data near '+eventname)
    # plt.savefig(eventname+'_L1_spectrogram.'+plottype)

if make_plots:
    #  plot the whitened data, zooming in on the signal region:

    # pick a shorter FTT time interval, like 1/16 of a second:
    NFFT = int(fs/16.0)
    # and with a lot of overlap, to resolve short-time features:
    NOVL = int(NFFT*15/16.0)
    # choose a window that minimizes "spectral leakage" 
    # (https://en.wikipedia.org/wiki/Spectral_leakage)
    window = np.blackman(NFFT)

    # Plot the H1 whitened spectrogram around the signal
    plt.figure(figsize=(10,6))
    spec_H1, freqs, bins, im = plt.specgram(strain_H1_whiten[indxt], NFFT=NFFT, Fs=fs, window=window, 
                                            noverlap=NOVL, cmap=spec_cmap, xextent=[-deltat,deltat])
    plt.xlabel('time (s) since '+str(tevent))
    plt.ylabel('Frequency (Hz)')
    plt.colorbar()
    plt.axis([-0.5, 0.5, 0, 500])
    plt.title('aLIGO H1 strain data near '+eventname)
    # plt.savefig(eventname+'_H1_spectrogram_whitened.'+plottype)

    # Plot the L1 whitened spectrogram around the signal
    plt.figure(figsize=(10,6))
    spec_H1, freqs, bins, im = plt.specgram(strain_L1_whiten[indxt], NFFT=NFFT, Fs=fs, window=window, 
                                            noverlap=NOVL, cmap=spec_cmap, xextent=[-deltat,deltat])
    plt.xlabel('time (s) since '+str(tevent))
    plt.ylabel('Frequency (Hz)')
    plt.colorbar()
    plt.axis([-0.5, 0.5, 0, 500])
    plt.title('aLIGO L1 strain data near '+eventname)
    # plt.savefig(eventname+'_L1_spectrogram_whitened.'+plottype)

# read in the template (plus and cross) and parameters for the theoretical waveform
try:
    f_template = h5py.File('LOSC_Event_tutorial/'+fn_template, "r")
except:
    print("Cannot find template file!")
    print("You can download it from https://losc.ligo.org/s/events/"+eventname+'/'+fn_template)
    print("Quitting.")
    # quit()  # kernel restart!
    
# extract metadata from the template file:
template_p, template_c = f_template["template"][...]
t_m1 = f_template["/meta"].attrs['m1']
t_m2 = f_template["/meta"].attrs['m2']
t_a1 = f_template["/meta"].attrs['a1']
t_a2 = f_template["/meta"].attrs['a2']
t_approx = f_template["/meta"].attrs['approx']
f_template.close()
# the template extends to roughly 16s, zero-padded to the 32s data length. The merger will be roughly 16s in.
template_offset = 16.

# whiten the templates:
template_p_whiten = whiten(template_p,psd_H1,dt)
template_c_whiten = whiten(template_c,psd_H1,dt)
template_p_whitenbp = filtfilt(bb, ab, template_p_whiten) / normalization
template_c_whitenbp = filtfilt(bb, ab, template_c_whiten) / normalization

# Compute, print and plot some properties of the template:

# constants:
clight = 2.99792458e8                # m/s
G = 6.67259e-11                      # m^3/kg/s^2 
MSol = 1.989e30                      # kg

# template parameters: masses in units of MSol:
t_mtot = t_m1+t_m2
# final BH mass is typically 95% of the total initial mass:
t_mfin = t_mtot*0.95

# Final BH radius, in km:
# 计算最终黑洞半径。
# 这里 G 是万有引力常数，t_mfin 是最终黑洞的质量（以太阳质量为单位），MSol 是太阳的质量，clight 是光速。
# 这个公式基于广义相对论中的施瓦茨希尔德半径（Schwarzschild radius），用于估算非旋转黑洞的半径。
# 物理上，它代表了黑洞的“事件视界”大小，即从该半径内部逃逸所需的逃逸速度大于光速。
R_fin = 2*G*t_mfin*MSol/clight**2/1000.

# complex template:
template = (template_p + template_c*1.j) 
ttime = time-time[0]-template_offset

# compute the instantaneous frequency of this chirp-like signal:
# 计算引力波的瞬时频率。
# tphase 是信号的相位，fs 是采样频率。
# 物理上，引力波的频率随着两个黑洞的合并过程而变化，这个频率提供了关于合并过程的重要信息。
# 瞬时频率的计算有助于分析引力波信号的特性，如频率随时间的变化。
tphase = np.unwrap(np.angle(template))
fGW = np.gradient(tphase)*fs/(2.*np.pi)
# fix discontinuities at the very end:
# iffix = np.where(np.abs(np.gradient(fGW)) > 100.)[0]
iffix = np.where(np.abs(template) < np.abs(template).max()*0.001)[0]
fGW[iffix] = fGW[iffix[0]-1]
fGW[np.where(fGW < 1.)] = fGW[iffix[0]-1]

# compute v/c:
# 引力波源的轨道速度与光速的比值。
# 这里 t_mtot 是系统的总质量，G 是万有引力常数，MSol 是太阳的质量，fGW 是引力波的频率，clight 是光速。
# 这个公式基于开普勒第三定律的修正形式，用于估算两个黑洞在其轨道上移动的速度。
# 在物理上，这个速度比值提供了关于黑洞合并过程中的动态信息，如系统的能量和动量变化。
voverc = (G*t_mtot*MSol*np.pi*fGW/clight**3)**(1./3.)

# index where f_GW is in-band:
f_inband = fband[0]
iband = np.where(fGW > f_inband)[0][0]
# index at the peak of the waveform:
ipeak = np.argmax(np.abs(template))

# number of cycles between inband and peak:
Ncycles = (tphase[ipeak]-tphase[iband])/(2.*np.pi)

print('Properties of waveform template in {0}'.format(fn_template))
print("Waveform family = {0}".format(t_approx))
print("Masses = {0:.2f}, {1:.2f} Msun".format(t_m1,t_m2))
print('Mtot = {0:.2f} Msun, mfinal = {1:.2f} Msun '.format(t_mtot,t_mfin))
print("Spins = {0:.2f}, {1:.2f}".format(t_a1,t_a2))
print('Freq at inband, peak = {0:.2f}, {1:.2f} Hz'.format(fGW[iband],fGW[ipeak]))
print('Time at inband, peak = {0:.2f}, {1:.2f} s'.format(ttime[iband],ttime[ipeak]))
print('Duration (s) inband-peak = {0:.2f} s'.format(ttime[ipeak]-ttime[iband]))
print('N_cycles inband-peak = {0:.0f}'.format(Ncycles))
print('v/c at peak = {0:.2f}'.format(voverc[ipeak]))
print('Radius of final BH = {0:.0f} km'.format(R_fin))

Properties of waveform template in GW150914_4_template.hdf5
Waveform family = lalsim.SEOBNRv2
Masses = 41.74, 29.24 Msun
Mtot = 70.98 Msun, mfinal = 67.43 Msun 
Spins = 0.35, -0.77
Freq at inband, peak = 43.05, 169.84 Hz
Time at inband, peak = -0.08, -0.02 s
Duration (s) inband-peak = 0.06 s
N_cycles inband-peak = 4
v/c at peak = 0.57
Radius of final BH = 199 km

if make_plots:
    plt.figure(figsize=(10,16))
    plt.subplot(4,1,1)
    plt.plot(ttime,template_p)
    plt.xlim([-template_offset,1.])
    plt.grid()
    plt.xlabel('time (s)')
    plt.ylabel('strain')
    plt.title(eventname+' template at D_eff = 1 Mpc')
    
    plt.subplot(4,1,2)
    plt.plot(ttime,template_p)
    plt.xlim([-1.1,0.1])
    plt.grid()
    plt.xlabel('time (s)')
    plt.ylabel('strain')
    
    plt.subplot(4,1,3)
    plt.plot(ttime,fGW)
    plt.xlim([-1.1,0.1])
    plt.grid()
    plt.xlabel('time (s)')
    plt.ylabel('f_GW')
    
    plt.subplot(4,1,4)
    plt.plot(ttime,voverc)
    plt.xlim([-1.1,0.1])
    plt.grid()
    plt.xlabel('time (s)')
    plt.ylabel('v/c')
    # plt.savefig(eventname+'_template.'+plottype)

from  tqdm import tqdm
def correlation_scratch(strain, template):
    result = np.zeros_like(strain)
    for i in tqdm(range(len(strain))):
        result[i] = np.dot(np.roll(np.append(template, np.zeros(template.size//2)),
                                   -template.size//2+i)[:strain.size], strain)
    return result

win = np.hanning(strain_H1_whitenbp.size)
# 数据
strain_scratch = strain_H1_whitenbp * win

# 模板
template_scratch = template_p_whitenbp * win

strain_scratch.shape, template_scratch.shape

((131072,), (131072,))

corr_scratch = correlation_scratch(strain_scratch, template_scratch)
corr_scratch.shape

100%|██████████| 131072/131072 [01:04<00:00, 2038.50it/s]

(131072,)

plt.subplot(3,1,1)
plt.plot(strain_scratch)
plt.subplot(3,1,2)
plt.plot(template_scratch)
plt.subplot(3,1,3)
plt.plot(abs(corr_scratch))

[<matplotlib.lines.Line2D at 0x7fdcd77cccd0>]

plt.subplot(3,1,1)
plt.plot(strain_scratch)
plt.xlim(6.5e4,6.8e4)
plt.subplot(3,1,2)
plt.plot(template_scratch)
plt.xlim(6.5e4,6.8e4)
plt.subplot(3,1,3)
plt.plot(abs(corr_scratch))
plt.xlim(6.5e4,6.8e4)

(65000.0, 68000.0)

index_slic_argmax = np.argmax(abs(corr_scratch))
index_slic_argmax

67340

i = index_slic_argmax

plt.figure(figsize=(8,4))
plt.plot(strain_scratch)
plt.ylim(-10,10)
template_scratch_rolled = np.roll(np.append(template_scratch, np.zeros(template_scratch.size//2)),
                                     -template_scratch.size//2+i)[:strain_scratch.size]
plt.twinx()
plt.plot(template_scratch_rolled, 'r')
plt.xlim(6.5e4,6.8e4)

(65000.0, 68000.0)

# Check
np.abs(np.dot(template_scratch_rolled , strain_scratch))

2035182.2367850973

# Check
max(abs(corr_scratch))

2035182.2367850973

# Numpy
# Cross-correlation of two 1-dimensional sequences.
corr_numpy = np.correlate(strain_scratch, template_scratch, mode='same')

# Check
np.allclose(corr_scratch, corr_numpy)

True

# Scipy
# Cross-correlate two N-dimensional arrays.
corr_scipy = scipy.signal.correlate(strain_scratch, template_scratch, mode='same')

# Check
np.allclose(corr_scratch, corr_scipy)

True

%timeit np.correlate(strain_scratch, template_scratch, mode='same')
%timeit scipy.signal.correlate(strain_scratch, template_scratch, mode='same')

1.83 s ± 9.74 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
10.7 ms ± 5.79 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# Numpy
# Cross-correlation of two 1-dimensional sequences.
conv_numpy = np.convolve(strain_scratch, template_scratch[::-1], mode='same')

# Check
np.allclose(corr_scratch, conv_numpy)

True

# Scipy
# Cross-correlate two N-dimensional arrays.
conv_scipy = scipy.signal.convolve(strain_scratch, template_scratch[::-1], mode='same')

# Check
np.allclose(corr_scratch, conv_scipy)

True

%timeit np.convolve(strain_scratch, template_scratch[::-1], mode='same')
%timeit scipy.signal.convolve(strain_scratch, template_scratch[::-1], mode='same')

1.82 s ± 10.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
11 ms ± 3.55 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

corr_fft = np.abs(np.fft.ifft( np.fft.fft(strain_scratch) * np.fft.fft(template_scratch).conjugate() )) 
print(corr_fft.shape)

plt.plot(corr_fft)

(131072,)

[<matplotlib.lines.Line2D at 0x7fdccaf7da80>]

# Check
corr_scipy = abs(scipy.signal.correlate(strain_scratch, template_scratch, mode='same'))
plt.plot(corr_scipy)
max(abs(corr_scipy)), np.argmax(abs(corr_scipy))

(2035182.2367850975, 67340)

# Roll the corr_fft
plt.plot(np.roll(corr_fft, strain_scratch.size//2))
max(corr_fft), np.argmax(np.roll(corr_fft, strain_scratch.size//2))

(2035182.237237697, 67340)

# Check
i = np.argmax(np.roll(corr_fft, strain_scratch.size//2))

plt.figure(figsize=(8,4))
plt.plot(strain_scratch)
plt.ylim(-10,10)
template_scratch_rolled = np.roll(np.append(template_scratch, np.zeros(template_scratch.size//2)),
                                     -template_scratch.size//2+i)[:strain_scratch.size]
plt.axvline(i, color='r', ls = '--')
plt.twinx()
plt.plot(template_scratch_rolled, 'r')
plt.xlim(6.5e4,6.8e4)

(65000.0, 68000.0)

image = plt.imread('LOSC_Event_tutorial/mf.png')
plt.figure(dpi=300)
plt.imshow(image)
plt.xticks([])
plt.yticks([]);

# -- To calculate the PSD of the data, choose an overlap and a window (common to all detectors)
#   that minimizes "spectral leakage" https://en.wikipedia.org/wiki/Spectral_leakage
NFFT = 4*fs
psd_window = np.blackman(NFFT)
# and a 50% overlap:
NOVL = NFFT//2

# define the complex template, common to both detectors:
template = (template_p + template_c*1.j) 
# We will record the time where the data match the END of the template.
etime = time+template_offset
# the length and sampling rate of the template MUST match that of the data.
datafreq = np.fft.fftfreq(template.size)*fs
df = np.abs(datafreq[1] - datafreq[0])

# to remove effects at the beginning and end of the data stretch, window the data
# https://en.wikipedia.org/wiki/Window_function#Tukey_window
try:   dwindow = scipy.signal.windows.tukey(template.size, alpha=1./8)  # Tukey window preferred, but requires recent scipy version 
except: dwindow = scipy.signal.windows.blackman(template.size)          # Blackman window OK if Tukey is not available

# prepare the template fft.
template_fft = np.fft.fft(template*dwindow) / fs

# loop over the detectors
dets = ['H1', 'L1']
for det in dets:

    if det == 'L1': data = strain_L1.copy()
    else:           data = strain_H1.copy()

    # -- Calculate the PSD of the data.  Also use an overlap, and window:
    data_psd, freqs = mlab.psd(data, Fs = fs, NFFT = NFFT, window=psd_window, noverlap=NOVL)

    # Take the Fourier Transform (FFT) of the data and the template (with dwindow)
    data_fft = np.fft.fft(data*dwindow) / fs

    # -- Interpolate to get the PSD values at the needed frequencies
    power_vec = np.interp(np.abs(datafreq), freqs, data_psd)

    # -- Calculate the matched filter output in the time domain:
    # Multiply the Fourier Space template and data, and divide by the noise power in each frequency bin.
    # Taking the Inverse Fourier Transform (IFFT) of the filter output puts it back in the time domain,
    # so the result will be plotted as a function of time off-set between the template and the data:
    optimal = data_fft * template_fft.conjugate() / power_vec
    optimal_time = 2*np.fft.ifft(optimal)*fs

    # -- Normalize the matched filter output:
    # Normalize the matched filter output so that we expect a value of 1 at times of just noise.
    # Then, the peak of the matched filter output will tell us the signal-to-noise ratio (SNR) of the signal.
    sigmasq = 1*(template_fft * template_fft.conjugate() / power_vec).sum() * df
    sigma = np.sqrt(np.abs(sigmasq))
    SNR_complex = optimal_time/sigma

    # shift the SNR vector by the template length so that the peak is at the END of the template
    peaksample = int(data.size / 2)  # location of peak in the template
    SNR_complex = np.roll(SNR_complex,peaksample)
    SNR = abs(SNR_complex)

    # find the time and SNR value at maximum:
    indmax = np.argmax(SNR)
    timemax = time[indmax]
    SNRmax = SNR[indmax]

    # Calculate the "effective distance" (see FINDCHIRP paper for definition)
    # d_eff = (8. / SNRmax)*D_thresh
    d_eff = sigma / SNRmax
    # -- Calculate optimal horizon distnace
    horizon = sigma/8

    # Extract time offset and phase at peak
    phase = np.angle(SNR_complex[indmax])
    offset = (indmax-peaksample)

    # apply time offset, phase, and d_eff to template 
    template_phaseshifted = np.real(template*np.exp(1j*phase))    # phase shift the template
    template_rolled = np.roll(template_phaseshifted,offset) / d_eff  # Apply time offset and scale amplitude
    
    # Whiten and band-pass the template for plotting
    template_whitened = whiten(template_rolled,interp1d(freqs, data_psd),dt)  # whiten the template
    template_match = filtfilt(bb, ab, template_whitened) / normalization # Band-pass the template
    
    print('For detector {0}, maximum at {1:.4f} with SNR = {2:.1f}, D_eff = {3:.2f}, horizon = {4:0.1f} Mpc' 
          .format(det,timemax,SNRmax,d_eff,horizon))

    if make_plots:

        # plotting changes for the detectors:
        if det == 'L1': 
            pcolor='g'
            strain_whitenbp = strain_L1_whitenbp
            template_L1 = template_match.copy()
        else:
            pcolor='r'
            strain_whitenbp = strain_H1_whitenbp
            template_H1 = template_match.copy()

        # -- Plot the result
        plt.figure(figsize=(10,8))
        plt.subplot(2,1,1)
        plt.plot(time-timemax, SNR, pcolor,label=det+' SNR(t)')
        #plt.ylim([0,25.])
        plt.grid('on')
        plt.ylabel('SNR')
        plt.xlabel('Time since {0:.4f}'.format(timemax))
        plt.legend(loc='upper left')
        plt.title(det+' matched filter SNR around event')

        # zoom in
        plt.subplot(2,1,2)
        plt.plot(time-timemax, SNR, pcolor,label=det+' SNR(t)')
        plt.grid('on')
        plt.ylabel('SNR')
        plt.xlim([-0.15,0.05])
        #plt.xlim([-0.3,+0.3])
        plt.grid('on')
        plt.xlabel('Time since {0:.4f}'.format(timemax))
        plt.legend(loc='upper left')
        # plt.savefig(eventname+"_"+det+"_SNR."+plottype)

        plt.figure(figsize=(10,8))
        plt.subplot(2,1,1)
        plt.plot(time-tevent,strain_whitenbp,pcolor,label=det+' whitened h(t)')
        plt.plot(time-tevent,template_match,'k',label='Template(t)')
        plt.ylim([-10,10])
        plt.xlim([-0.15,0.05])
        plt.grid('on')
        plt.xlabel('Time since {0:.4f}'.format(timemax))
        plt.ylabel('whitened strain (units of noise stdev)')
        plt.legend(loc='upper left')
        plt.title(det+' whitened data around event')

        plt.subplot(2,1,2)
        plt.plot(time-tevent,strain_whitenbp-template_match,pcolor,label=det+' resid')
        plt.ylim([-10,10])
        plt.xlim([-0.15,0.05])
        plt.grid('on')
        plt.xlabel('Time since {0:.4f}'.format(timemax))
        plt.ylabel('whitened strain (units of noise stdev)')
        plt.legend(loc='upper left')
        plt.title(det+' Residual whitened data after subtracting template around event')
        # plt.savefig(eventname+"_"+det+"_matchtime."+plottype)
                 
        # -- Display PSD and template
        # must multiply by sqrt(f) to plot template fft on top of ASD:
        plt.figure(figsize=(10,6))
        template_f = np.absolute(template_fft)*np.sqrt(np.abs(datafreq)) / d_eff
        plt.loglog(datafreq, template_f, 'k', label='template(f)*sqrt(f)')
        plt.loglog(freqs, np.sqrt(data_psd),pcolor, label=det+' ASD')
        plt.xlim(20, fs/2)
        plt.ylim(1e-24, 1e-20)
        plt.grid()
        plt.xlabel('frequency (Hz)')
        plt.ylabel('strain noise ASD (strain/rtHz), template h(f)*rt(f)')
        plt.legend(loc='upper left')
        plt.title(det+' ASD and template around event')
        # plt.savefig(eventname+"_"+det+"_matchfreq."+plottype)

For detector H1, maximum at 1126259462.4395 with SNR = 18.6, D_eff = 814.44, horizon = 1889.6 Mpc
For detector L1, maximum at 1126259462.4324 with SNR = 13.2, D_eff = 999.74, horizon = 1650.6 Mpc

Amplitude spectra	Fourier transform	Power spectra
Two-sided (peak)	$\frac{A}{2} = \frac{A_{RMS}}{\sqrt 2}$	$\frac{A^2}{4} = \frac{{A_{RMS}}^2}{2}$
Two-sided (RMS)	$\frac{A}{2 \cdot \sqrt 2} = \frac{A_{RMS}}{2}$	$\frac{A^2}{8} = \frac{{A_{RMS}}^2}{4}$
One-sided (peak)	$A = \sqrt 2 \cdot A_{RMS}$	$\frac{A^2}{2} = {A_{RMS}}^2$
One-sided (RMS)	$\frac{A}{\sqrt 2} = A_{RMS}$	$\frac{A^2}{4} = \frac{{A_{RMS}}^2}{2}$

Binary Black Hole Signals in LIGO Open Data¶

Intro to signal processing¶

Download the data on a computer with a python installation¶

Set the event name to choose event and the plot type¶

Read the event properties from a local json file (download in advance at `LOSA_Event_tutorial/`):¶

Read in the data¶

Data Gaps¶

First look at the data from H1 and L1¶

【Basic】Understanding the FFT Algorithm by Python¶

采样定理与 Nyquist 频率¶

傅里叶变换¶

FFT from scratch (`DFT_slow`)¶

FFT from scratch (`FFT`)¶

FFT from scratch (`FFT_vectorized`)¶

Fast FFT from scratch (FFT_fast)¶

Verification of correctness for the code¶

Efficiency of the algorithms¶

频域索引和频域坐标之间的映射关系 (`numpy.fft.fftfreq`)¶

Amplitude Spectrum¶

Power Spectrum¶

Windowing¶

【Basic】Spectral Analysis (谱分析)¶

PSD¶

Specgram¶

Plot the Amplitude Spectral Density (ASD)¶

Whitening + bandpass¶

Spectrograms¶

Waveform Template¶

不得不略过的知识点¶

我就是想看看匹配度，整那么多数学干啥？ (Naive 匹配滤波)¶

互相关（Cross-correlation）¶

卷积（Convolution）¶

卷积定理（Convolution theorem）¶

Matched filtering to find the signal¶

Binary Black Hole Signals in LIGO Open Data¶

Intro to signal processing¶

Download the data on a computer with a python installation¶

Set the event name to choose event and the plot type¶

Read the event properties from a local json file (download in advance at LOSA_Event_tutorial/):¶

Read in the data¶

Data Gaps¶

First look at the data from H1 and L1¶

【Basic】Understanding the FFT Algorithm by Python¶

采样定理与 Nyquist 频率¶

傅里叶变换¶

FFT from scratch (DFT_slow)¶

FFT from scratch (FFT)¶

FFT from scratch (FFT_vectorized)¶

Fast FFT from scratch (FFT_fast)¶

Verification of correctness for the code¶

Efficiency of the algorithms¶

频域索引和频域坐标之间的映射关系 (numpy.fft.fftfreq)¶

Amplitude Spectrum¶

Power Spectrum¶

Windowing¶

【Basic】Spectral Analysis (谱分析)¶

PSD¶

Specgram¶

Plot the Amplitude Spectral Density (ASD)¶

Whitening + bandpass¶

Spectrograms¶

Waveform Template¶

不得不略过的知识点¶

我就是想看看匹配度，整那么多数学干啥？ (Naive 匹配滤波)¶

互相关（Cross-correlation）¶

卷积（Convolution）¶

卷积定理（Convolution theorem）¶

Matched filtering to find the signal¶

Read the event properties from a local json file (download in advance at `LOSA_Event_tutorial/`):¶

FFT from scratch (`DFT_slow`)¶

FFT from scratch (`FFT`)¶

FFT from scratch (`FFT_vectorized`)¶

频域索引和频域坐标之间的映射关系 (`numpy.fft.fftfreq`)¶